Cod diet: spatial analysis of cod and flounder diets in relation to density
Max Lindmark & Michele Casini
12 November, 2021
Intro
Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod
# Load libraries, install if needed
library(tidyverse); theme_set(theme_classic(base_size = 10))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
library(bayesplot)
library(tidybayes)
library(brms)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB) # To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/cod_flounder_diets_spatial_cache/html")For maps
# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
ggplot(swe_coast_proj) + geom_sf() 
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Make default base map plot
plot_map_raster <-
ggplot(swe_coast_proj) +
geom_sf(size = 0.3) +
labs(x = "Longitude", y = "Latitude") +
theme_facet_map(base_size = 14)
# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}
theme_clean <- function() {
theme_minimal(base_family = "Barlow Semi Condensed") +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(family = "BarlowSemiCondensed-Bold"),
axis.title = element_text(family = "BarlowSemiCondensed-Medium"),
strip.text = element_text(family = "BarlowSemiCondensed-Bold",
size = rel(1), hjust = 0),
strip.background = element_rect(fill = "grey80", color = NA))
}Read and plot data
# These data are for total and prey specific stomach models
cod <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> time_period = col_character(),
#> quarter_fact = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
cod <- cod %>%
mutate(year = as.integer(year),
quarter = as.factor(quarter),
depth2_sc = depth - mean(depth),
saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
filter(year > 2014) %>%
filter(!quarter == 2) %>%
drop_na(predfle_density_sc, predcod_density_sc) %>%
droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#> converted 'quarter' from double to factor (0 new NA)
#> new variable 'depth2_sc' (double) with 106 unique values and 0% NA
#> new variable 'saduria_entomon_per_mass' (double) with 1,709 unique values and 0% NA
#> new variable 'tot_prey_biom_per_mass' (double) with 7,202 unique values and 0% NA
#> new variable 'depth_sc' (double) with 106 unique values and 0% NA
#> filter: removed 4,968 rows (57%), 3,703 rows remaining
#> filter: removed 91 rows (2%), 3,612 rows remaining
#> drop_na: no rows removed
fle <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/fle_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
fle <- fle %>%
mutate(year = as.integer(year),
quarter = as.factor(quarter),
depth2_sc = depth - mean(depth),
saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
filter(!quarter == 2) %>%
drop_na(predfle_density_sc, predcod_density_sc) %>%
droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#> converted 'quarter' from double to factor (0 new NA)
#> new variable 'depth2_sc' (double) with 81 unique values and 0% NA
#> new variable 'saduria_entomon_per_mass' (double) with 608 unique values and 0% NA
#> new variable 'tot_prey_biom_per_mass' (double) with 1,735 unique values and 0% NA
#> new variable 'depth_sc' (double) with 81 unique values and 0% NA
#> filter: removed 58 rows (3%), 2,190 rows remaining
#> drop_na: no rows removed
# These data are for the diet analysis (all prey groups)
cod_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/pca_cod_diet_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "COD") %>% filter(year > 2014) %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 4,968 rows (57%), 3,703 rows remaining
#> filter: removed 91 rows (2%), 3,612 rows remaining
fle_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/pca_fle_diet_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "FLE") %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 58 rows (3%), 2,190 rows remainingRead the prediction grids
# And now read in pred_grid2 which has oxygen values at location and time and depth:
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid2.csv") ### CHANGE TO OUR NEW PRED GRID
# Clean
pred_grid2 <- pred_grid2 %>%
mutate(year = as.integer(year)) %>%
drop_na(depth)
# Add ices_rect
pred_grid2$ices_rect <- ices.rect2(pred_grid2$lon, pred_grid2$lat)
# pred_grid2_q1 <- pred_grid2 %>% mutate(quarter = factor(1))
# pred_grid2_q4 <- pred_grid2 %>% mutate(quarter = factor(4))Calculate Schoener’s overlap index for diet
# Plot data in space
plot_map_raster +
geom_point(data = fle_prey, aes(x = X * 1000, y = Y * 1000, color = "FLE"), alpha = 0.5) +
geom_point(data = cod_prey, aes(x = X * 1000, y = Y * 1000, color = "COD"), alpha = 0.5) +
theme_plot()
cod_prey <- cod_prey %>% filter(lat < 58)
#> filter: removed 19 rows (1%), 3,593 rows remaining
fle_prey <- fle_prey %>% filter(lat < 58)
#> filter: removed 10 rows (<1%), 2,180 rows remaining
# Reformat data to calculate Schoeners overlap index
# colnames(fle_prey)
fle_prey_long <- fle_prey %>%
pivot_longer(14:29) %>%
group_by(name, year, quarter, ices_rect) %>% # ices rect?
summarise(fle_diet = sum(value)) %>%
arrange(name, year, ices_rect) %>%
ungroup() %>%
group_by(year, quarter, ices_rect) %>%
mutate(fle_diet_tot = sum(fle_diet),
fle_diet_prop = fle_diet / fle_diet_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x30, now 34880x16]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,296 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> mutate (grouped): new variable 'fle_diet_tot' (double) with 81 unique values and 0% NA
#> new variable 'fle_diet_prop' (double) with 445 unique values and 1% NA
#> ungroup: no grouping variables
fle_prey_wide <- fle_prey_long %>%
pivot_wider(names_from = name, values_from = fle_diet_prop, values_fill = 0)
#> pivot_wider: reorganized (name, fle_diet_prop) into (amphipoda_tot, bivalvia_tot, clupea_harengus_tot, clupeidae_tot, gadiformes_tot, …) [was 1296x7, now 526x21]
cod_prey_long <- cod_prey %>%
pivot_longer(14:29) %>%
group_by(name, year, quarter, ices_rect) %>%
summarise(cod_diet = sum(value)) %>%
arrange(name, year, ices_rect) %>%
ungroup() %>%
group_by(year, quarter, ices_rect) %>%
mutate(cod_diet_tot = sum(cod_diet),
cod_diet_prop = cod_diet / cod_diet_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x30, now 57488x16]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,696 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> mutate (grouped): new variable 'cod_diet_tot' (double) with 106 unique values and 0% NA
#> new variable 'cod_diet_prop' (double) with 784 unique values and 0% NA
#> ungroup: no grouping variables
cod_prey_wide <- cod_prey_long %>%
pivot_wider(names_from = name, values_from = cod_diet_prop, values_fill = 0)
#> pivot_wider: reorganized (name, cod_diet_prop) into (amphipoda_tot, bivalvia_tot, clupea_harengus_tot, clupeidae_tot, gadiformes_tot, …) [was 1696x7, now 892x21]
# Calculate Schoener index
schoener <- left_join(cod_prey_long, fle_prey_long) %>%
drop_na(name) %>%
drop_na(fle_diet_prop) %>%
drop_na(cod_diet_prop) %>%
group_by(year, quarter, ices_rect) %>%
summarise(schoener = 1 - 0.5*(sum(abs(fle_diet_prop - cod_diet_prop)))) %>%
ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect")
#> left_join: added 3 columns (fle_diet, fle_diet_tot, fle_diet_prop)
#> > rows only in x 432
#> > rows only in y ( 32)
#> > matched rows 1,264
#> > =======
#> > rows total 1,696
#> drop_na: no rows removed
#> drop_na: removed 448 rows (26%), 1,248 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 78 rows and 4 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables
# Calculate Levin niche index
levin <- left_join(cod_prey_long, fle_prey_long) %>%
drop_na(name) %>%
drop_na(fle_diet_prop) %>%
drop_na(cod_diet_prop) %>%
group_by(year, quarter, ices_rect) %>%
summarise(levin_cod = (1/(16-1)) * (((1/sum(cod_diet_prop^2))) - 1),
levin_fle = (1/(16-1)) * (((1/sum(fle_diet_prop^2))) - 1)) %>%
ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect")
#> left_join: added 3 columns (fle_diet, fle_diet_tot, fle_diet_prop)
#> > rows only in x 432
#> > rows only in y ( 32)
#> > matched rows 1,264
#> > =======
#> > rows total 1,696
#> drop_na: no rows removed
#> drop_na: removed 448 rows (26%), 1,248 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 78 rows and 5 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables
# Merge the indicies
ind <- left_join(levin, schoener)
#> Joining, by = c("year", "quarter", "ices_rect")
#> left_join: added one column (schoener)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 78
#> > ====
#> > rows total 78
# Summarise cod and flounder data by ices_rect then add to diet data
colnames(cod)
#> [1] "FR_tot" "FR_sad"
#> [3] "FR_spr" "saduria_entomon_tot"
#> [5] "tot_prey_biom" "common_tot"
#> [7] "unique_pred_id" "year"
#> [9] "quarter" "time_period"
#> [11] "quarter_fact" "pred_weight_g"
#> [13] "pred_cm" "predator_spec"
#> [15] "predcod_density" "predfle_density"
#> [17] "predcod_density_sc" "predfle_density_sc"
#> [19] "depth" "X"
#> [21] "Y" "lat"
#> [23] "long" "ices_rect"
#> [25] "cruise" "depth2_sc"
#> [27] "saduria_entomon_per_mass" "tot_prey_biom_per_mass"
#> [29] "depth_sc"
cod$year_rect_id <- paste(cod$year, cod$quarter, cod$ices_rect, sep = "_")
dens_sum <- cod %>% group_by(year_rect_id) %>%
summarise(predfle_density_tot = sum(predfle_density),
predcod_density_tot = sum(predcod_density)) %>%
ungroup() %>%
mutate(predfle_density_tot_sc = (predfle_density_tot - mean(predfle_density_tot)) / sd(predfle_density_tot),
predcod_density_tot_sc = (predcod_density_tot - mean(predcod_density_tot)) / sd(predcod_density_tot))
#> group_by: one grouping variable (year_rect_id)
#> summarise: now 112 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'predfle_density_tot_sc' (double) with 110 unique values and 0% NA
#> new variable 'predcod_density_tot_sc' (double) with 110 unique values and 0% NA
ind$year_rect_id <- paste(ind$year, ind$quarter, ind$ices_rect, sep = "_")
ind <- left_join(ind, dens_sum)
#> Joining, by = "year_rect_id"
#> left_join: added 4 columns (predfle_density_tot, predcod_density_tot, predfle_density_tot_sc, predcod_density_tot_sc)
#> > rows only in x 0
#> > rows only in y (34)
#> > matched rows 78
#> > ====
#> > rows total 78
# Summarise depth from the prediction grid then add to diet data
pred_grid2_sum <- pred_grid2 %>%
group_by(ices_rect) %>%
summarise(mean_depth = mean(depth)) %>%
mutate(depth_sc = (mean_depth - mean(mean_depth)) / sd(mean_depth)) %>%
ungroup()
#> group_by: one grouping variable (ices_rect)
#> summarise: now 61 rows and 2 columns, ungrouped
#> mutate: new variable 'depth_sc' (double) with 61 unique values and 0% NA
#> ungroup: no grouping variables
ind <- left_join(ind, dplyr::select(pred_grid2_sum, ices_rect, depth_sc, mean_depth))
#> Joining, by = "ices_rect"
#> left_join: added 2 columns (depth_sc, mean_depth)
#> > rows only in x 0
#> > rows only in y (43)
#> > matched rows 78
#> > ====
#> > rows total 78
# Plot the indicies
ggplot(ind) +
geom_jitter(aes(factor(year), schoener),
alpha = 0.8, width = 0.2, height = 0, size = 2)
ggplot(ind) +
geom_jitter(aes(factor(year), levin_cod, color = "cod"),
alpha = 0.8, width = 0.2, height = 0, size = 2) +
geom_jitter(aes(factor(year), levin_fle, color = "fle"),
alpha = 0.8, width = 0.2, height = 0, size = 2) +
scale_color_brewer(palette = "Dark2")
# Against response variables
# Schoener
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = schoener)) +
geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predfle_density_tot, y = schoener)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predcod_density_tot, y = schoener)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
# Now Levin
# Cod
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = levin_cod)) +
geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 
ggplot(ind, aes(x = predfle_density_tot, y = levin_cod)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 
ggplot(ind, aes(x = predcod_density_tot, y = levin_cod)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 
# Flounder
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = levin_fle)) +
geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 
ggplot(ind, aes(x = predfle_density_tot, y = levin_fle)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 
ggplot(ind, aes(x = predcod_density_tot, y = levin_fle)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter) 
# Add refs here or in text already! THis is complementary to prey weights in stomach
# Do cumsum plots
# Now split by size-classes, check Orior co author paper for which sizes to split
# Think about which model to fit to these indices?
# ---- brms with random rectangle and fixed year? Evaluate differences in intercept and plot covariates?
# ... Then when I have size-group? I get like 5 "species". Maybe do pca on those?
# Do that after resolving old comments, adding new diet data and correctly adding quarter and new cpue data!!!
# Quickly check data to determine which distribution to use
ind %>%
ungroup() %>%
count(schoener == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#> `schoener == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 75 0.962
#> 2 TRUE 3 0.0385
ind %>%
ungroup() %>%
count(levin_cod == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#> `levin_cod == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 77 0.987
#> 2 TRUE 1 0.0128
ind %>%
ungroup() %>%
count(levin_fle == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#> `levin_fle == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 75 0.962
#> 2 TRUE 3 0.0385
# Fit beta models, so few zeroes!
ind %>% arrange(schoener) %>% dplyr::select(schoener)
#> # A tibble: 78 × 1
#> schoener
#> <dbl>
#> 1 0
#> 2 0
#> 3 0
#> 4 0.000490
#> 5 0.00127
#> 6 0.00134
#> 7 0.00149
#> 8 0.00209
#> 9 0.00240
#> 10 0.00261
#> # … with 68 more rows
ind %>% arrange(levin_cod) %>% dplyr::select(levin_cod)
#> # A tibble: 78 × 1
#> levin_cod
#> <dbl>
#> 1 0
#> 2 0.00173
#> 3 0.00265
#> 4 0.00493
#> 5 0.00801
#> 6 0.0132
#> 7 0.0163
#> 8 0.0179
#> 9 0.0191
#> 10 0.0270
#> # … with 68 more rows
ind %>% arrange(levin_fle) %>% dplyr::select(levin_fle)
#> # A tibble: 78 × 1
#> levin_fle
#> <dbl>
#> 1 0
#> 2 0
#> 3 0
#> 4 0.000321
#> 5 0.000737
#> 6 0.00191
#> 7 0.00206
#> 8 0.00304
#> 9 0.00779
#> 10 0.00811
#> # … with 68 more rows
ind <- ind %>% mutate(schoener2 = ifelse(schoener == 0, 0.0001, schoener),
levin_cod2 = ifelse(levin_cod == 0, 0.0001, levin_cod),
levin_fle2 = ifelse(levin_fle == 0, 0.0001, levin_fle),
year_f = as.factor(year),
quarter_f = as.factor(quarter),
ices_rect_f = as.factor(ices_rect))
#> mutate: new variable 'schoener2' (double) with 76 unique values and 0% NA
#> new variable 'levin_cod2' (double) with 78 unique values and 0% NA
#> new variable 'levin_fle2' (double) with 76 unique values and 0% NA
#> new variable 'year_f' (factor) with 6 unique values and 0% NA
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> new variable 'ices_rect_f' (factor) with 18 unique values and 0% NAFit brms models with densities as covariates to diversity and overlap indices
# All covariates
m_schoen_beta_full <- brm(
bf(schoener2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f + quarter_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Trying to compile a simple C file
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Users/maxlindmark/Library/R/4.0/library/Rcpp/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/unsupported" -I"/Users/maxlindmark/Library/R/4.0/library/BH/include" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#> ^
#> ;
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include <complex>
#> ^~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1
#> Start sampling
plot(m_schoen_beta_full)
conditional_effects(m_schoen_beta_full)
loo_m_schoen_beta_full <- loo(m_schoen_beta_full, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta_full)
# Only density covariates
m_schoen_beta_dens <- brm(
bf(schoener2 ~ 0 + year_f + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Start sampling
plot(m_schoen_beta_dens)
conditional_effects(m_schoen_beta_dens)
loo_m_schoen_beta_dens <- loo(m_schoen_beta_dens, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta_dens)
# Simplest model
m_schoen_beta <- brm(
bf(schoener2 ~ 0 + year_f + (1|ices_rect_f),
phi ~ 0 + year_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Start sampling
plot(m_schoen_beta)
conditional_effects(m_schoen_beta)
loo_m_schoen_beta <- loo(m_schoen_beta, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta)
loo_compare(loo_m_schoen_beta, loo_m_schoen_beta_dens, loo_m_schoen_beta_full)
#> elpd_diff se_diff
#> m_schoen_beta 0.0 0.0
#> m_schoen_beta_dens -1.0 1.4
#> m_schoen_beta_full -4.3 2.6
# Plot models (working with the m_schoen_beta_dens)
summary(m_schoen_beta_dens)
#> Family: beta
#> Links: mu = logit; phi = log
#> Formula: schoener2 ~ 0 + year_f + predfle_density_tot_sc + predcod_density_tot_sc + (1 | ices_rect_f)
#> phi ~ 0 + year_f
#> Data: ind (Number of observations: 78)
#> Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
#> total post-warmup draws = 12000
#>
#> Group-Level Effects:
#> ~ices_rect_f (Number of levels: 18)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.19 0.16 0.01 0.58 1.00 4846 5380
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> year_f2015 -2.74 0.46 -3.59 -1.75 1.00 10133
#> year_f2016 -2.91 0.30 -3.48 -2.29 1.00 11525
#> year_f2017 -2.57 0.36 -3.24 -1.83 1.00 11899
#> year_f2018 -2.14 0.51 -3.08 -1.08 1.00 10237
#> year_f2019 -1.58 0.51 -2.55 -0.55 1.00 12042
#> year_f2020 -1.46 0.42 -2.26 -0.60 1.00 11745
#> predfle_density_tot_sc -0.25 0.16 -0.57 0.06 1.00 10524
#> predcod_density_tot_sc 0.18 0.15 -0.11 0.48 1.00 10331
#> phi_year_f2015 2.36 0.64 0.95 3.46 1.00 9842
#> phi_year_f2016 2.30 0.38 1.49 2.97 1.00 11640
#> phi_year_f2017 1.82 0.42 0.91 2.56 1.00 12310
#> phi_year_f2018 1.21 0.54 0.02 2.17 1.00 10486
#> phi_year_f2019 0.67 0.51 -0.40 1.58 1.00 11405
#> phi_year_f2020 1.46 0.50 0.40 2.35 1.00 12856
#> Tail_ESS
#> year_f2015 7731
#> year_f2016 7531
#> year_f2017 8457
#> year_f2018 7513
#> year_f2019 8082
#> year_f2020 8680
#> predfle_density_tot_sc 9258
#> predcod_density_tot_sc 8973
#> phi_year_f2015 7032
#> phi_year_f2016 7508
#> phi_year_f2017 7677
#> phi_year_f2018 7340
#> phi_year_f2019 8729
#> phi_year_f2020 8334
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
# Working with the posterior
posterior_beta <- m_schoen_beta_dens %>%
gather_draws(`b_.*`, regex = TRUE) %>%
mutate(component = ifelse(str_detect(.variable, "phi_"), "Precision", "Mean"),
intercept = str_detect(.variable, "Intercept"))
ggplot(posterior_beta, aes(x = .value, y = fct_rev(.variable), fill = component)) +
geom_vline(xintercept = 0) +
stat_halfeye(aes(slab_alpha = intercept), alpha = 0.5,
.width = c(0.8, 0.95), point_interval = "median_hdi") +
scale_fill_brewer(palette = "Dark2") +
scale_slab_alpha_discrete(range = c(1, 0.4)) +
guides(fill = "none", slab_alpha = "none") +
labs(x = "Coefficient", y = "Variable") +
facet_wrap(vars(component), ncol = 1, scales = "free_y") +
theme_plot() +
NULL
# Selected parameters
# https://www.sciencedirect.com/topics/mathematics/logit-scale
posterior_beta %>%
filter(component == "Mean" & .variable %in% c("b_predcod_density_tot_sc", "b_predfle_density_tot_sc")) %>%
ggplot(., aes(x = plogis(.value), y = .variable, fill = .variable)) +
geom_vline(xintercept = 0) +
stat_halfeye(alpha = 0.5, .width = c(0.8, 0.95)) +
scale_fill_brewer(palette = "Dark2") +
labs(x = "Coefficient", y = "Variable") +
theme_plot() +
guides(fill = FALSE) +
NULL
#> filter (grouped): removed 144,000 rows (86%), 24,000 rows remaining
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
# Prediction! (non dens model)
beta_bayes_pred_yr <- m_schoen_beta %>%
epred_draws(newdata = tibble(year_f = as.factor(c(2015:2020))),
re_formula = NA)
ggplot(beta_bayes_pred_yr, aes(x = .epred, y = year_f, color = year_f, fill = year_f)) +
stat_halfeye(.width = c(0.8, 0.95), alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
guides(fill = FALSE, color = FALSE) +
coord_flip(xlim = c(0, 0.4)) +
labs(x = "Predicted Schoener's overlap index", y = NULL) +
theme_plot()
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
beta_bayes_pred_re <- m_schoen_beta %>%
epred_draws(newdata = tibble(expand.grid(year_f = as.factor(c(2015:2020)),
ices_rect_f = as.factor(unique(ind$ices_rect)))))
ggplot(beta_bayes_pred_re, aes(x = .epred, y = ices_rect_f, color = ices_rect_f, fill = ices_rect_f)) +
stat_halfeye(.width = c(0.8, 0.95), alpha = 0.7) +
scale_fill_viridis(discrete = TRUE) +
scale_color_viridis(discrete = TRUE) +
guides(fill = FALSE, color = FALSE) +
coord_cartesian(xlim = c(0, 0.4)) +
#coord_flip(xlim = c(0, 0.4)) +
labs(x = "Predicted Schoener's overlap index", y = NULL) +
theme_plot()
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
# Cod Levin index: all covariates
m_levin_beta_full_cod <- brm(
bf(levin_cod2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f + quarter_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 2, iter = 2000, warmup = 1000, cores = 2, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_levin_beta_full_cod)
conditional_effects(m_levin_beta_full_cod)
loo_m_schoen_beta_full <- loo(m_levin_beta_full_cod, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(m_levin_beta_full_cod)
m_levin_beta_full_fle <- brm(
bf(levin_fle2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f + quarter_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 2, iter = 2000, warmup = 1000, cores = 2, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_levin_beta_full_fle)
